Mapping directional variations in seismic character using gray-level co-occurrence matrix-based attributes

Texture attributes describe the spatial arrangement of neighboring amplitudes values within a given analysis window. We chose a statistical texture classification method, the gray-level co-occurrence matrix (GLCM), and its derived attributes, to produce a semiautomated description of the spatial arrangement of seismic facies. The GLCM is a measure of how often different combinations of neighboring pixel values occur. We tested the application of directional GLCM-based attributes for the detection of seismic variability within paleoriver features. Calculation of 3D GLCM-based attributes can be done in 13 space directions. The results of GLCM-based attribute calculation differed depending on the chosen GLCM parameters (number of gray levels, analysis window, and direction of calculation). We specifically focused on how the direction of calculation influenced the computation of attributes, while keeping other parameters constant. We first tested the workflow on a 2D training image and later ran on a real seismic amplitude volume from the Vienna Basin. Based on the GLCM-based attributes, we could map the channel features and extract them as geobodies. Additionally, we generated a new set of directional GLCM-based attributes to detect spatial changes in the seismic facies. By comparing these directional attributes, we could determine areas within the channel features having higher directional variability. Areas with higher tendency to directional variations might be associated with changes in lithology, seismic facies, or with seismic anisotropy. Introduction Common seismic attributes such as complex trace (Taner et al., 1979), coherence (Bahorich and Farmer, 1995), curvature (Marfurt and Kirlin, 2000), or spectral decomposition attributes (Partyka et al., 1999) use mathematical formulations to capture the geometry or physical properties of the subsurface and can be used to illuminate geologic features of interest. Texture is defined by the spatial configuration of rock units and is more diagnostic of and relevant to deformational fabrics, depositional facies, and reservoir properties than an averaged acoustic property (Gao, 2011). The methodology follows the way a seismic interpreter analyzes seismic amplitudes and waveforms. Among many methods available for texture analysis, we choose a statistical texture classification method known as the gray-level co-occurrence matrix (GLCM) (Haralick et al., 1973) and its derived attributes to produce a semiautomated description of the spatial arrangement of seismic amplitude values. The GLCM was primarily designed for texture classification of 2D images. It is a method widely used in many application domains, such as image classification of satellite imagery (e.g., Franklin et al., 2001; Tsai et al., 2007), sea-ice images (e.g., Soh and Tsatsoulis, 1999; Maillard et al., 2005), magnetic resonance, or computed tomography images (e.g., Kovalev et al., 2001; Zizzari et al., 2011). This methodology has been applied to seismic data for the past 20 years (Vinther et al., 1996; Gao, 1999, 2003, 2007, 2008a, 2008b, 2009, 2011; West et al., 2002; Chopra and Alexeev, 2005, 2006a, 2006b; Yenugu et al., 2010; de Matos et al., 2011; Eichkitz et al., 2013). To calculate GLCMbased attributes on 3D seismic data, it is necessary to adapt the method to work in a 3D space. The number of possible directions increases in comparison with the 2D images. For this research project, we used four horizontal directions (0°, 45°, 90°, and 135°) to determine the directional behavior of seismic amplitude values. The workflow was tested on a data set from the Vienna Basin, Austria; the objective was to illuminate seismic character variability within channel features. Method The GLCM is a measure of how often different combinations of neighboring pixel values occur in an image. For a 2D image, the immediate neighboring pixels can Joanneum Research, Geophysics and Geothermics, Leoben, Austria. E-mail: [email protected]; marcellus.schreilechner@ joanneum.at; [email protected]. dGB Earth Sciences, Enschede, Netherlands. E-mail: [email protected]. Manuscript received by the Editor 14 May 2014; revised manuscript received 30 July 2014; published online 30 December 2014. This paper appears in Interpretation, Vol. 3, No. 1 (February 2015); p. T13–T23, 9 FIGS. http://dx.doi.org/10.1190/INT-2014-0099.1. © 2014 Society of Exploration Geophysicists and American Association of Petroleum Geologists. All rights reserved. t Technical papers Interpretation / February 2015 T13 D ow nl oa de d 01 /1 2/ 15 to 1 43 .2 24 .2 53 .2 34 . R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / be analyzed in four different directions (0°, 45°, 90°, and 135°). The four directions can also be combined to form an average GLCM. Such computation can eliminate the influences of bed dip and azimuth to a certain degree (Gao, 2007). In 3D data, the number of possible directions increases to thirteen. To illustrate this case, the Rubik’s cube is used (Figure 1a–1c), which is built of 27 smaller cubes. The cube in the center, the turning point in a Rubik’s Cube, is the reference point for which the calculations are performed (Figure 1b). The center point is surrounded by 26 neighboring cubes. If we take the center point and draw lines from it to all neighboring cubes, we get 13 directions (see Figure 1c). As in the 2D case, it is possible to calculate the GLCM in single directions, combine several directions, or calculate an average GLCM. Additionally, our algorithm integrates apparent local structural dip of the data (dip steering, e.g., de Rooij and Tingdahl, 2002). With dip steering, the GLCM input volume is warped along the seismic stratigraphy by following a precalculated 3D dip field. One advantage of dip steering is that the input volume for the GLCM computation does not mix signals from different seismic packages, leading to sharper images that are easier to interpret. For further information on the most common methods to estimate structural dip, we recommend the original works by Barnes (1996), Luo et al. (1996), Marfurt et al. (1998), Bakker (2002), Hoecker and Fehmers (2002), and Tingdahl and de Groot (2003). The first step in GLCM calculation is the transformation of the seismic amplitude volume (usually stored in 32 bit floating point) into a gray-level cube. For seismic data, 16 (4 bit) or 32 (5 bit) gray levels are usually regarded as sufficient. Calculations with a greater number of gray levels do not result in any significant differences in the computed attributes (Chopra and Alexeev, 2006a; Gao, 2007). Calculations with more gray levels result in only minor improvements in resulting images and are costly because of long computation times. Despite these issues, we believe that a higher number of gray levels are important to improve the signal-to-noise ratio (S/N), offering a better delineation of the seismic facies. Enhancement of the S/N can be observed up to approximately 512 (9 bit) gray levels. After that point, only minor changes are observable. We integrate a linked-list approach (Clausi and Jernigan, 1998; Clausi and Zhao, 2002; Clausi and Zhao, 2003) to overcome computational performance issues. A linked list is a sequence of nodes in which each node consists of a data part and a reference to the next node. Using this approach, nonzero matrix entries are skipped and calculation time is no longer a function of the number of gray levels. Only the size of the analysis window affects computation times. From the GLCM, a series of textural attributes can be derived. Haralick et al. (1973) originally introduce 14 attributes. Since then, a few more attributes have been developed and published (e.g., Soh and Tsatsoulis, 1999; Wang et al., 2010). Prior to computing any GLCM-based attributes, the GLCM is transformed by Figure 1. (a) In a 3D case, the number of neighbors for one sample point can be best explained by examining a Rubik’s Cube. (b) The center of the Rubik’s Cube (the red box in panel [b]) has in total 26 neighboring boxes. The boxes are aligned in 13 directions. (c) Analogous to this, a sample point in a seismic volume has 26 neighbors aligned in 13 directions. In the developed workflow, it is possible to calculate GLCM-based attributes along a single direction, a combination of directions (e.g., inline direction, crossline direction, etc.), or all directions at once (after Eichkitz et al., 2013). T14 Interpretation / February 2015 D ow nl oa de d 01 /1 2/ 15 to 1 43 .2 24 .2 53 .2 34 . R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / dividing each entry by the sum of all entries, expressing it as a kind of probability matrix. The GLCM-based attributes can be divided into three general groups. First, the contrast group includes attributes such as contrast, homogeneity, and cluster tendency (Haralick et al., 1973; Wang et al., 2010). The attributes in this group are basically a function of the probability of each matrix entry and the difference between gray levels (i and j). Therefore, these contrast group attributes are related to the distance from the GLCM diagonal. Values close to the diagonal (where i and j are the same) result in lower contrast, whereas the contrast increases as the distance from the diagonal increases. In the sample image (Figure 2), contrast is lower for the 135° direction (Figure 2g) than in the 45° direction (Figure 2f). For the 135° direction, the matrix entries are aligned around the diagonal where i and j are the same. Second, the orderliness group includes attributes such as energy and entropy (Haralick et al., 1973). Attributes in the orderliness group measure how regularly gray-level values are distributed within a given searchwindow. Unlike the contrast group, attributes in the orderliness group are solely a function of the GLCM probability entries. Third, the statistics group includes attributes such as mean and variance (Haralick et al., 1973). These are the common mean and variance calculated from the GLCM probabilities. We generated a small 2D sample image to demonstrate the GLCM calculation along different directions (Figure 2a). In this example, a northwest–southeast-trending feature was modeled. Gray levels along the main feature Figure 2. (a) Example for the calculation of GLCM-based attributes using eight gray levels from a synthetic 2D gray-scale image. (b) The gray-scales of the image correspond to discrete values. (c) The number of co-occurrences of pixel pairs for a given search window is counted and a GLCM is generated. In this synthetic case, the analysis window is 7 × 7 and attribute response would be assigned to the center image point. Based on this co-occurrence matrix, several attributes can be calculated. In this example, the GLCM are determined for the (d) 0°, (e) 90°, (f) 45°, and (g) 135° directions and (h) for all directions at once. The first step in calculation is the determination of co-occurrences (figure parts in the middle). Zero entries are marked in gray, and the highest value of each matrix is marked in dark green. It is evident that calculations in single directions lead to sparse matrices. The GLCM is normalized by the sum of the entries to get a kind of probability matrix (figures on the right). Finally, the probabilities are used for the calculation of the GLCM-based attributes. In the column labeled “attribute response,” the results for entropy, contrast, homogeneity, entropy, and cluster tendency are shown. Interpretation / February 2015 T15 D ow nl oa de d 01 /1 2/ 15 to 1 43 .2 24 .2 53 .2 34 . R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / are similar, and the largest variation occurs perpendicular to the main feature. Eight gray levels (3 bit) were used (Figure 2b). The number of gray levels (Ng) defines the rank of the square matrix (Figure 2c). In the example, the GLCM is a square matrix of rank eight. In Figure 2d–2g, the results for GLCM calculations for four directions are shown. Additionally, the results of the combinedGLCMcalculation in all four possible directions are shown in Figure 2h. Note the high number of zero values inside the matrix, which heavily influences computation times.With the use of a linked-list approach, we can overcomethisproblem.ThecalculatedGLCM-basedattributes energy, homogeneity, and cluster tendency are highest when the calculation is done parallel (135° direction) to the main feature in the tested image. Contrast gives the highest values perpendicular to themain feature. Entropy is smallest in the strike (45°) direction, whereas in other directions, entropy values are higher but similar to each other. From these results, we can see that the attribute response is dependent on the strike of the main feature. The directional behavior of the GLCM-based attributes can thus be used to estimate variability of the seismic amplitude responses. The input to generate GLCM-based attributes can be seismic amplitude volumes, or any derived seismic attribute volume (coherence, curvature, spectral decomposition). Against this background, we developed a workflow composed of four steps (Figure 3). Currently, the workflow for mapping directional variations operates in 4 directions only, but we intend to extend this to a full analysis of all 13 possible space directions. The first step is the calculation of GLCM-based attributes in four separate space directions (0°, 45°, 90°, and 135°). By integrating dip steering methodology into GLCM calculation, we can increase the S/N. In the second step, we determine for each sample point and each GLCM-based attribute the maximum and minimum values. Additionally, we store the direction in which maximum and minimum values occur for each attribute. Information of maximum and minimum values and their directions are stored in separate cubes. Subsequently, we use this information to identify the direction of highFigure 3. The workflow for seismic amplitude variability detection by using GLCMbased attributes is divided into four parts. (a) The first part is the calculation of each GLCM-based attribute in four space directions (in this example, only the four horizontal directions are used). (b) The second step is the determination of minima and maxima values and their direction for each GLCM-based attribute. (c) In step 3, the ratio between maximum and minimum is calculated. The ratio is used to determine a threshold value. (d) The threshold value is applied in step 4 to identify areas with higher directional variability and to visualize only these directions of minima and maxima. T16 Interpretation / February 2015 D ow nl oa de d 01 /1 2/ 15 to 1 43 .2 24 .2 53 .2 34 . R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / est variation in the seismic data. Rarely do we get the same attribute response in all directions, although minimum and maximum values sometimes are very close to each other. In these areas, we would overestimate the directional variability in our seismic data. Therefore, it is necessary to condition the output data to distinguish between areas with high-directional variability and areas with low directional variability. This is done in steps three and four of the workflow. In step 3, we calculate the ratio between maximum and minimum values and set a manually determined threshold value for each attribute. The threshold value is then used in step 4 to set all areas with a ratio below the threshold value to have no directional variability. In areas above the threshold, the minimum and maximum directions are determined. By combining the calculated GLCM-based attributes, we are able to determine areas that tend to have lower directional variability from those areas prone to have more directional variability. Application to a channel system within the Vienna Basin We apply the workflow to a 3D seismic data set from the Vienna Basin (close to the border between Austria and Slovakia). The bin size of the cube is 25 × 25 m (82 × 82 ft). At a depth of approximately 1000 ms, a meandering channel system was interpreted. The channel deposits are most probably of Sarmatian or Pannonian age (Middle to Upper Miocene), and have a width of 130–300 m (426–984 ft) and an average thickness of 70 ms. Strauss et al. (2006) describe similar features in the southern Vienna Basin. The channel system is clearly identified with the help of the coherence attribute (Eichkitz et al., 2012). We interpret the edges of the channel as lineaments on several time slices (Figure 4) and use the interpreted lineaments as guidance for the interpretation of the GLCM-based attribute displays. For comparison, the GLCM-based attributes (energy, homogeneity, contrast, entropy, and homogeneity) were computed in four directions (0°, 45°, 90°, and 135°) with and without dip steering applied to the input seismic amplitude volume. Sixty-four gray levels were used for GLCM-based calculation. In the energy and cluster tendency attributes, the channel system is well imaged (Figure 5) because the surrounding rocks show very low, uniform values and the channel filling has high values in both attributes. In the homogeneity and entropy attributes, the channel system is still visible, but it is more difficult to differentiate the channel from the surrounding strata. The homogeneity attribute is more difficult to interpret because the surrounding areas have many spots with responses similar to that inside the channel system. For the contrast attribute, the interpretation of the channel system is very difficult. No real differentiation between the channel and the surrounding rocks is possible. As noted, all attributes were calculated in four horizontal space directions with and without dip steering. In general, the results using dip steering have a higher S/N (Figure 6), so we decided to use the results with dip steering for the determination of spatial variation in seismic data. In the energy attribute display, large variations are observed in all four space directions. Visual Figure 4. In the tested data set, multiple channel features are observed. To visualize these channel features and to interpret the channel edges, coherence time slices were used (Eichkitz et al., 2012). In panel (a) is a time slice of the amplitude and in panel (b) is a time slice of semblance-based coherence. The interpretation of the channel edges is displayed in panel (c). Coherence interpretation is overlain on all GLCM-based attribute plots. Interpretation / February 2015 T17 D ow nl oa de d 01 /1 2/ 15 to 1 43 .2 24 .2 53 .2 34 . R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / analysis of these four different images (Figure 7) leads to the conclusion that calculation in the inline direction gives the highest values for this attribute. According to theory and supported by the test image (Figure 2), highenergy values can be interpreted as the direction in which the least variation occurs. The 135° direction (northwest–southwest) generally shows the lowest energy response. This can be interpreted as the direction in which the highest variation in seismic data response is present. Ratio cubes for each attribute volume were generated in the next step. For each GLCM-based attribute, the maximum value is divided by the minimum value. Based on these ratio cubes, we manually set the threshold values. For energy and entropy, a threshold value of 1.2 is best for illuminating the channel interior. For Figure 5. GLCM-based attributes calculated along 0° direction (east–west): (a) Energy, (b) contrast, (c) homogeneity, (d) entropy, and (e) cluster tendency. The channel structure can clearly be seen in the energy and cluster tendency attribute. The contrast attribute in this direction reveals little information about the channel structure. In panel (f) , the coherence time slice is shown once more for comparison reasons. T18 Interpretation / February 2015 D ow nl oa de d 01 /1 2/ 15 to 1 43 .2 24 .2 53 .2 34 . R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg / homogeneity and cluster tendency, a threshold value of 1.1 was used. Using these values, we produce the final images of direction of maximum and minimum GLCMbased attribute values. Figure 8 shows the direction of maximum values. Within the channel features, the maximum values mostly occur in the 0° and 90° direcFigure 7. GLCM-based energy attribute calculated in four different space directions. In panel (a), the zero direction (east–west), in (b), the 90° direction (north–south), in (c), the 45° direction (northeast–southwest), and in (d), the 135° direction (northwest–southeast) are shown. Figure 6. An important step in seismic attribute calculation is to take the structural dip into account. For this purpose, we determine the structural dip by using a complex trace analysis, a discrete scan, or a gradient structure tensor approach. In this work, the discrete scan was found to produce best results. In panels (a and c), the images represent the calculation without using the structural dip. In panels (a and b), several seismic traces are shown. Clearly, one can see that the vertical positions of peak amplitudes are not aligned in one single horizontal line. By neglecting the structural dip of the data, the center of our analysis window would fall at the dark gray line in panel (a) and our analysis window would be a perfect cube. Calculation of GLCM-based energy in the 0° direction (c) without steering and (d) with steering. In the steered computation, the S/N is higher and certain features can better be imaged (see the gray arrows). Interpretation / February 2015 T19 D ow nl oa de d 01 /1 2/ 15 to 1 43 .2 24 .2 53 .2 34 . R ed is tr ib ut io n su bj ec t t o SE G li ce ns e or c op yr ig ht ; s ee T er m s of U se a t h ttp :// lib ra ry .s eg .o rg /